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^ ■ Abstract 

' Coupled backward and forward wave amplitudes of an electromagnetic field propagating in a periodic 

and nonlinear medium at Bragg resonance are governed by the nonlinear coupled mode equations (NL- 
CME). This system of PDEs, similar in structure to the Dirac equations, has gap soliton solutions that 
travel at any speed between and the speed of light. A recently considered strategy for spatial trapping 
or capture of gap optical soliton light pulses is based on the appropriate design of localized defects in the 
, periodic structure. Localized defects in the periodic structure give rise to defect modes, which persist as 

nonlinear defect modes as the amplitude is increased. Soliton trapping is the transfer of incoming soliton 
energy to nonlinear defect modes. To serve as targets for such energy transfer, nonlinear defect modes 
must be stable. We therefore investigate the stability of nonlinear defect modes. Resonance among 
discrete localized modes and radiation modes plays a role in the mechanism for stability and instability, 
in a manner analogous to the nonlinear Schrodinger / Gross-Pitaevskii (NLS/GP) equation. However, 
the nature of instabilities and how energy is exchanged among modes is considerably more complicated 
' than for NLS/GP due, in part, to a continuous spectrum of radiation modes which is unbounded above 

If) I and below. In this paper we (a) establish the instability of branches of nonlinear defect states which, for 

vanishing amplitude, have a linearization with eigenvalue embedded within the continuous spectrum, (b) 
numerically compute, using Evans function, the linearized spectrum of nonlinear defect states of an in- 
00 ' teresting multiparameter family of defects, and (c) perform direct time-dependent numerical simulations 

, in which we observe the exchange of energy among discrete and continuum modes. 
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1 Introduction 

In optical fiber communication systems, data is carried as pulses of light. Expensive and rate-limiting 
steps in these systems come in processing the data at so-called optical / electrical interfaces. Processing 
is typically done electronically; signals are converted to electronic form, read, processed, and retransmit- 
ted. A major goal, therefore, is to bypass the optical / electrical interface and implement "all-optical" 
processing by using the nonlinear optical properties of the medium. Hence, there has been great interest 
in finding novel materials and optical structures (arrangements of materials) which eff'ect light in different 
ways. Among these arc Bragg gratings and Bragg gratings with defects-one dimensional arrangements, 
as well as higher-dimensional structures such as photonic crytal fibers in which the grating structure is 
transverse to the direction of propagation [43], and other two and three-dimensional structures [16]. 

In an optical fiber Bragg grating, the refractive index of the glass varies periodically, with period 
resonant with the carrier wavelength of the propagating light; see figure 1.1. Light propagation through 
such fibers has several interesting properties that makes them useful in technological applications. The 
grating structure couples forward-moving light at the resonant wavelength to backward-moving light of 
the same wavelength. In the low-amplitude (linear) limit, this makes the fiber opaque to light in a certain 
range of wavelengths, the so-called photonic band-gap. When the amplitude of the light is increased, 
the band-gap is shifted due to the nonlinear dependence of refractive index on intensity. Thus, a range 
of wavelengths that are non-propagating at low intensity are shifted into the range of propagating 
wavelengths (the pass-band) at higher intensities. This is the mechanism by which localized pulses 
known as gap solitons exist; see subsection 2.2. In the regime of weak nonlincarity and Bragg resonance, 
Maxwell's equations can be reduced, via multiple scale asymptotic methods, to the nonlinear coupled 
mode equations (NLCME), reviewed in section 2.1 [2, 15, 19, 32].^. 

Gap solitons may, in theory, propagate (in the stationary reference frame) at any speed between 
and c/n. Here, c denotes the vacuum speed of light, n the refractive index of the optical fiber core and 
c/n is the speed of light in the fiber without the grating. In [29], we showed via modeling, analysis and 
numerical simulation, that propagating optical gap solitons, could be trapped at specially-constructed 
defects in the grating structure. Trapping and interaction of gap solitons in a number of related structures 
is considered in [46, 47, 14] . Recent experimental advances have made possible the slowing of propagating 
gap soliton fight pulses from 0.5 x c [23] to 0.16 x c [52]. 

""^A special case of NLCME (equation (2.1) with the self-phase modulation terms omitted) is the massive Thirring system [25], 
a completely integrable PDE modeling the interaction of massive fermions. The methods described in this article apply equally 
well to the massive Thirring system. 
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Figure 1.1: Schematic of Bragg resonance condition. Periodic refractive index (lower figure) and electric 
field envelope with carrier wave in Bragg resonance (A = 2 x period). 

The focus of the present paper is on the stabihty and dynamics of such trapped hght. Locahzed 
defects in a grating appear as spatially localized potentials in the NLCME model. In the linear (low 
light intensity) limit, the resulting linear coupled mode equations with potentials have spatially localized 
linear defect eigenstates; see subsection 2.3. As the intensity is increased from zero, nonlinear defect 
modes, "pinned" at the defect location, bifurcate from the zero state at the linear eigenfrequencies; see 
subsection 2.4. 

Trapping of a gap soliton by a defect can be understood as the resonant transfer of energy from 
an incoming soliton to a pinned nonlinear defect mode. In [29] we demonstrated through numerical 
experiments such resonant energy transfer / trapping, for sufficiently slow soliton pulses. The analogous 
question has also been studied for the nonlinear Schrodinger / Gross-Pitaevskii (NLS/GP) equations; 
see [28, 35, 36, 34] and references cited therein. For the purpose of comparison, we review the stability 
of and interactions between nonlinear defect modes of NLS/GP in section 3. 

In order for the energy localized in a defect to remain spatially confined in a nonlinear defect mode, 
it is necessary that the mode be stable. Thus, in this paper we consider the stability of nonlinear defect 
modes for a large class of defects introduced in [29] by 

• studying the linearized spectral problem about different families (branches in the global bifurcation 
diagram) of nonlinear defect modes; see sections 4 and 5 for analytical perturbation theory and 

numerics, and 

• studying time-dependent simulations of the initial value problem for NLCME. We consider de- 
fects which support multiple nonlinear linear defect modes. As the time-evolution proceeds these 
nonlinear bound state families compete for the energy localized in the defect; see section 6. 

The linearized stability of a nonlinear defect mode is governed by a spectral problem of the form: 

= /3 V', = -S3, H* = H; (1.1) 

see section 4. Eigenvalues are complex values of P, for which (1.1) has a nontrivial solution, giving 
rise to a time-dependent solution of the linearized dynamics (j) — V'e~'^*. The self-adjoint operator, H, 
can be expressed as H = Hq + W , where Hq corresponds to the linear (zero intensity) coupled mode 
equation operator and W tends to zero quadratically in the amplitude (L^-norm) of the nonlinear defect 
mode about which we have linearized. 

The continuous spectrum of 'EsH typically consists of two symmetric semi-infinite real intervals with 
an open gap centered at /3 = 0. A priori, since is not self-adjoint, discrete eigenvalues may lie 

anywhere in the complex plane, subject to constraints inherited by (1.1) from NLCME, a Hamiltonian 
system. In particular, if /3 is an eigenvalue, then so are —(3,(3* and —(3*. Thus, a necessary condition 
for linearized stability of the flow idt<j> — ^sHcp is that the spectrum of Y.qH be real. 

Now Ti^Ho has stable spectrum and may have real discrete eigenvalues in the spectral gap or real 
embedded eigenvalues within the continuous spectrum. Numerous different instability-onset scenarios 
corresponding to different types of bifurcations arise as Sgi/o deforms to Ti^H = S3(i/o + W), as the 
intensity {L^ norm) of a nonlinear defect mode is increased along the bifurcation branch. We focus on 
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Figure 1.2: Instability onset scenario 1: as the amplitude of the nonlinear defect mode is increased from 
to e the embedded frequencies ib/^o split into complex frequencies ib/3g it i'^^. 



one such scenario here and in more detail in section 4.2. The other scenarios are discussed briefly in 
section 4.3. 

The first scenario occurs when Ssffo has a symmetric pair of real embedded eigenvalues within the 
continuous spectrum. Generically these perturb, for || W|| positive and arbitrary, to two pairs of complex 
conjugate eigenvalues and therefore yield instability. The mechanism is related to the notion of Krein 
signature. A detailed perturbation calculation demonstrating this instability is presented in section 4.2, 
showing that the instability is of order or equivalently, of order equal to the fourth power of 

the defect mode amplitude. In the particular examples we study numerically, these instability rates are 
observed, but are seen in some cases to be quite small. 

In section 6 we explore the temporal dynamics of the initial value problem for NLCME via direct 
numerical simulation. We present simulations with various defects supporting one. two and three modes 
which illustrate both the spectral instability scenarios explored in section 4 and the temporal dynamics 
of energy exchange among defect modes and radiation modes. The mechanisms are similar to, but 
more complicated than, those studied by SofFcr and Weinstein [67, 68] for the nonlinear Schrodingcr / 
Gross-Pitaevskii (NLS / GP) equation, who consider NLS / GP with a "defect potential" supporting 
two bound states and thus, two branches of nonlinear defect modes (ground and excited), which compete 
for energy confined by the potential; sec also [71]. In the NLS/GP problem, resonance of an embedded 
eigenvalue, associated with the linear excited bound state, with the continuous spectrum leads to energy 
transfer from the excited state to the ground state and to radiation modes. The asymptotic distribution 
of excited state energy is also studied; see [67, 68, 74]. 

We recap the structure of the paper. In section 2, we introduce NLCME, its gap solitons, and two 
families of defects and their linear and nonlinear bound states. We also summarize the results of a 
previous paper in which these defects arc used to trap gap solitons. In section 3 we review some anal- 
ogous results for the NLS equation with a localized defect. Section 4 discusses the properties of the 
linearized operator, outlining several types of instability with particular attention to the bifurcations of 
emhedded frequencies. Section 5 provides a brief introduction to the Evans function, a useful analytical 
tool for exploring the discrete spectrum of the linearized operator, as well numerical results using the 
Evans function to study the spectral stability of nonlinear defect modes. Section 6 shows time-dependent 
simulations confirming the stability predictions of the previous section. Section 7 contains a short sum- 
mary and discussion. The appendices contain a list of symbols, a derivation of NLCME with potentials 
from the problem of wave propagation in a Bragg grating with defects, and a detailed description the 
numerical measures that are taken to ensure accurate computations of both the nonlinear defect modes 
and the discrete spectrum of the linearization. 



2 NLCME and NLCME with defects 



2.1 The Nonlinear Coupled Mode Equations 

Consider an electromagnetic wave-packet, a pulse which is spectrally concentrated about a carrier wave- 
length, propagating in a waveguide whose refractive index varies periodically in the direction of prop- 
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agation. If carrier wavelength is in resonance with the waveguide periodicity (Bragg resonance), then 
backward and forward waves couple strongly. 

Envelope equations for these forward and backward wave amplitudes satisfy the nonlinear coupled 
mode equations. In appendix B we show that localized defects in the periodic structure give rise to the 
following modification of the NLCME, which we call by the same name, as derived in [29, 73]: 

idTE+ + idzE+ + k{Z)E_ + V{Z)E+ + (|£;+|2 + 2\E_\'')E+ = 
idj.E_ - idzE_ + K.{Z)E+ + V{Z)E_ + (|£;_|2 + 2\E+\^)E_ =0. 

The coefficient functions ("potentials") k and V are determined from the refractive index profile. Hero. 
Z is the "slow" coordinate along the direction of propagation, T is a slow time variable and the full 
electric field is given by 

E = E+{Z, T)e'(^-*) + i;_(Z,T)e-'(^+*) + c.c, 

i.e. E+ and i?_ arc complex envelopes of rapidly varying electromagnetic fields, assumed to be small 
amplitude and linearly polarized. The length and time variables [z, t) are chosen such that the wavenum- 
ber and frequency are both one, i.e. if in dimensional variables {z,i), the carrier wave has wavelength 
27r/A and period 27rc/(nA), then we choose z = 2ttz/X and t = 2'Kci/{nX). The condition that the 
grating is uniform (exactly periodic) away from the defect region implies that k{Z) Koo, constant, 
and V{Z) ^ as \Z\ oo. Defining E — {^'^), equation (2.1) may be rewritten as 

{idr + iaadz + V{Z) + k{Z)(jx) E+J\f{E, E*)E = 0, (2.2) 

Where ai and are the Pauli matrices ci = ( 5 J ) and 173= (0-1)1 the superscript asterisk represents 
complex conjugation, and TV represents the nonlinear term of (2.1), 

M[E,E)-y Q \E_\^ + 2\E+\y- ^^■'^> 

A detailed derivation, including careful accounting of the nondimensionalization and scalings, is given 
in [29, 32, 73]. 

The NLCME system (2.1) has two conserved quantity, the total intensity I{E) and the Hamiltonian 
(energy) H{E). The total intensity, I{E), is defined as 



I{E)^ {\E+\' + \E.\')dZ, (2.4) 
J —00 

i.e. the square of the norm of the solution. The Hamiltonian is given by 

H{E) = f (iE+dzE+ - iE_dzE_ + k{Z) {E+E*_ + E^EJ) 
J —00 

+ V{Z) {\E+\^ + \E_\^) + 2\E+\''\E_\^ + + dZ (2.5) 

In the absence of a defect {k{Z) constant, V{Z) = 0), NLCME also conserves momentum, as a conse- 
quence of Noether's theorem. 

The primary focus of this paper is the standing wave solutions of (2.2) of the form 

E{Z,T)=£{Z)e-''''^, 

for which (2.2) reduces to the nonlinear eigenvalue problem 

{oj + ia^dz + V{Z) + K{Z)ai) £ + Af{£, £*)£ = 0, £ £ x L^. (2.6) 
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2.2 Constant coefficient NLCME — The spectral gap and "gap solitons" 

NLCME with a uniform grating has k{Z) = k^, constant, and V{Z) = 0, and has been studied 
extensively. In the hnearized system, setting the cubic terms to zero, one may look for solutions of the 
form 



E.J' ' ' Ve- 

and find the dispersion relation = +A:^, so that for frequencies |a;| < k^o, k is purely imaginary, i.e. 
there exists a spectral gap (— Koo, Koo)- This results in the reflection of light at the resonant frequencies 
in the gap. 

The fully nonlinear coupled mode system was shown in [2, 15] to support uniformly propagating 
bound states called gap solitons, parameterized by a velocity v and a detuning parameter 6 satisfying 
\v\ < 1 and < 5 < tt: 

E± = ± A^^a^^ (sin 6) e'(''+") sech {0 T iS/2) ; (2.7) 

where: 



i=i — 1 ;„ = 



2(1 -t;2) / e2^ + e-'* 



z{smS){z — vt) ; a = —z^^—-{cos6){vz — t). 



In the Lorenz-shifted reference frame, the term e**^ is responsible for an internal oscillation with 
frequency k^o cos S which is inside the spectral gap, the nonlinearity having served to effectively shift 
the gap. This is also known as self-induced transparency. A comprehensive introduction to gap solitons 
is given in the review paper by de Sterke and Sipe [19]. Gap solitons may propagate in the laboratory 
reference frame at any velocity v below the speed of light (here c = 1 by the leading to B.6) and thus are 
intriguing as possible components of all-optical communications systems. Zero-speed waves, for example, 
are candidates for bits in an optical buffer or memory device. 

2.3 Defect modes of the Uneair coupled mode equations 

In the zero-intensity limit, we ignore the nonlinear term in (2.6), giving a linear eigenvalue problem 

( CO J + hv,n ) = i.'^J + lo^idz + V{Z) + K{Z)ai) £^ = 0, 

£*£L^ X L^. (2.8) 

Hereafter, frequencies and vectors with subscript asterisks will refer to solutions the linear eigenprob- 
lem (2.8), and those without will refer to solutions of (2.6). 

Note that hy.K, is self-adjoint on x and therefore has real spectrum. Moreover, we assume 
V{Z) — + and n{Z) Kqo rapidly as \Z\ oo, and therefore this spectral problem has two branches 
of continuous spectrum and a finite number of discrete frequencies given by 

spec {ujj + hv,^) = spec^„^ti^„„„3 U spec^ji^^^^^ = {w, G M : > Koo} U {wj* : j = 1, . . . , iV}, 

where the discrete frequencies satisfy — Kqo < < i^oo- 

We consider two families of defects in the current study which we refer to as "even defects" and "odd 
defects" . 

(a) Even defects 

If we specify k,{Z) to be even and V{Z) = 0, then it is simple to show that if a;» is a real eigenvalue 
of (2.8), then so is — w*. A simple argument eliminates the possiblity of nontrivial null eigenstates. Thus, 
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system (2.8) must have an even number of discrete eigenvalues with corresponding LF' eigenfunctions. 
As we do not know of any k{Z) for which the eigenvalue problem (2.8) may be solved in closed form, 
we choose a relatively simple form of k(-), 

k(Z) = 1 - 6sech(A:Z). (2.9) 

Thus, = 1 a-nd the linear operator oj^I + /iq.k has continuous spectrum {a; G M : |cij| > 1}. We have 
found via numerical simulation that as b is increased with k held fixed, the existing frequencies migrate 
toward the origin, and new discrete frequencies bifurcate in pairs from opposite endpoints of the band 
gap at uo = ±1. Such emergence from the edges has been studied for related problems, for example, 
in [42, 58]. 

(b) Odd defects 

In [29], we construct a three-parameter family ^ of defects of the form 

Lonk"^ sech^ [kZ) 
2{w'^ + n'^k'^ tanh^ (kZ)) " 



.{Z) = Vo;^ + n-k- tanh^ (kZ); V{Z) = .,.^.^,^\.' - (2-10) 



with fc > 0, n > 0, and w nonzero. ^ Standing wave solutions exist of the form: 

/ e'® \ 

^(0)* = ( _gg_i0 J sech" {kZ)e-"^°''^ , wo* = w, s = signw (2.11) 

/■^x./.x 1 nfctanh(fcZ) 
e = / V{0 dC= - arctan ^ ^. (2.12) 

For n > 1, the defect supports a total of 2[nJ — 1 distinct eigenvalues, where [nj is the greatest integer 
less than or equal to n. Its eigenvalues are 

LO±j. = ± + (2nj - f)k^; i = 1, . . . , n - 1. (2.13) 

If Lo and k are hold constant, while n is increased, the "ground state" eigenvalue loq* ~ ^ remains 
constant, while the other eigenvalues move from the edge of the spectral gap toward iwo*- Each time n 
passes through an integer value, two new eigenvalues are created in edge bifurcations from the two ends 
of the continuous spectrum. 

Reference [31] contains the general formula for the eigenfunctions. The solutions can ultimately 
be expressed in terms of hypcrgeometric functions of tanh kZ, which simplify to Legendre functions of 
tanh kZ when n G Z, and can be expressed as algebraic combinations of tanhfcZ and sechkZ. The 
general formula is quite complicated. 

In much of what follows we specialize to the case n = 2, in which case there are three linear defect 
modes. The first eigenfunction, with frequency uiq = w, is given by (2.11). The remaining two, with 
frequencies w±i* = ±-\/tJM-3fc^, are given (in a non-normalized form) by 

= ( ^t^f^'*^TT.f:ZL) sech^^. (2.14) 
' \3 {k - i{uj±i^ + Lo) tanh kZ) e-^'^(^) J ^ ' 

2.4 Defect modes of nonlinear coupled mode equation 

Consider equation (2.6), the nonlinear eigenvalue problem solved by the nonlnear defect modes. In 

the small amplitude (linear) limit, we expect solutions to be well-approximated by those of the linear 
eigenvalue problem (2.8), whose explicit eigenstates are displayed in section (2.3). 



^By a rescaling of the space and time coordinates, we may set A: = 1, so in reality this is a two-parameter family, but the 
above form makes it easier to create different defects with the same limiting grating strength Koo- 

^In the special case w = 0, then V{Z) = and k{Z) = nfctanh (kZ). The expressions for the eigenmodes are also altered. 
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a) 



Figure 2.1: The intensity as a function of the frequency of the numerically calculated nonlinear defect 
modes, with parameter values lo = 1, k = 1, n = 2. The band edges are at zt-y/S and the Hnear eigenvalues 
are chq = 1 and lj±i = ±2. The darker curve is plots the amplitude and frequency of a stationary gap 
soliton (2.7). 



In [29], wc used perturbation analysis to show that in analogy with the NLS/GP case [60], there 
exist nonlinear defect mode solutions of (2.6) 

= «(4>(Z)+0(|a|2)) and 
u = LJj^ + q\af + 0(|a|''), a e C, a ^ 

Notice that q < 0, due to the focusing character of the nonlinearity. Thus, as the amplitude is increased, 
the nonhnear frequency is shifted toward the left edge of the spectral gap. 

We aim to compute the nonlinear modes and their associated frequencies as accurately as possible. 
This is important so that errors in the numerical solutions contribute as little as possible to errors in 
their numerically calculated linearized spectrum and determination of stability. Details of the numerical 
calculation are provided in appendix C. Briefly, the derivatives are discretized using Fourier transforms, 
while the infinite interval is truncated using a nonuniform-grid method and the resulting algebraic 
equations derived are solved using MINPACK routines [53]. 

Figure 2.1 shows a typical example of the dependence of the defect mode's frequency on its amplitude. 
The frequency of each of the three defect modes shifts to the left as the intensity is increased; g < 0. 
At small intensities the frequency depends linearly on the intensity, and at larger intensities, the curves 
steepen near the left band edge. We believe that these branches terminate at the left band edge, but 
the nonlinear modes decay very slowly as the edge is approached, and thus become difficult to calculate 
numerically. 

2.5 Summciry of previous numerical experiments 

In [29], we explored the behavior of gap solitons incident on localized defects of the form (2.10). We 
found parameter regimes in which the soliton — or at least the energy it carries — can be captured. We 
hypothesized the existence of a nonlinear resonance mechanism between the soliton, and a nonlinear 
defect mode. By (2.7). the gap soliton oscillates with an internal frequency of about KoqCOsS as it 
propagates. If there exists a nonlinear defect mode of the same frequency and smaller L^-norm, then 
the solitary wave may transfer its energy to the defect mode. In particular, as 5 is increased from to 



(2.15) 
(2.16) 

(2.17) 
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Figure 2.2: The projections onto the U2 (thin soHd hne) and ui (dashed), showing that energy is originally, 
and quickly transferred from the moving soliton to the u)2 mode, and then slowly transferred to the uji 
mode. The thick line shows the total L/^ norm of the solution, which is not conserved due to absorbing 
boundary conditions at the endpoints. 



TT this internal frequency decreases from to — Koo- This is the thickest curve in figure 2.1; a family 
of states of the translation invariant NLCME, which bifurcates from the zero state at the band edge 
frequency. 

For sufficiently small 5 (frequency near the right band edge, marked (a) in the figure), there exists no 
nonlinear defect mode of the same frequency and lower amplitude to which the gap soliton can transfer 
its energy. In this case, the defect behaves as a barrier. For solitary waves above a critical velocity, 
the pulse passes by the defect with almost all its amplitude and its original velocity almost unchanged. 
Below the critical velocity, the solitary wave is reflected, again almost elastically. If, by contrast, there 
exists a nonlinear defect mode resonant and of with the solitary wave and of smaller amplitude, then 
the pulse may be trapped, for example points (b) and (d). In this case there exists a critical velocity, 
below which solitary waves are transmitted (inelastically) and below which they are largely trapped. 
Behavior at the point marked (c) is similar to that at (a)-the defect mode that bifurcates from = 2 
is a larger-amplitude state than the gap soliton, so it cannot be excited, and the defect mode bifurcation 
from Wo* is not resonant with this frequency. 

Similar behavior, namely the dichotomy between the nonrcsonant and elastic transmit/reflect behav- 
ior and the resonant and inelastic capture/transmit behavior, has also been seen for nonlinear Schrodinger 
solitons [28, 44]. Estimates for Vc in several related problems and an explanation for the existence of a 
critical velocity arc given in [27] and references tlica'cin. This trapping phenomenon was subsequently 
seen by Dohnal and Aceves for a two-dimensional generalization of equation (2.1) [1, 21]. 

Trapping could be more accurately described as a transfer of energy from the traveling soliton mode 
to a stationary nonlinear defect mode. This is pictured in figure 2.2 reprinted from [29]. This figure 
describes the evolution of a pulse captured by a wide defect supporting 5 linear defect modes. Energy 
initially moves into the mode associated with the frequency 0^2 but is subsequently transferred to the 
mode associated with frequency Wi. A part of the energy is transferred to radiation modes and is 
dispersed to infinity. It is thus important for us to understand the stability of these nonlinear modes, 
to assess their suitability as long-time containers for electromagnetic energy, as well as the mechanisms 
through which energy moves among discrete and continuum modes. 

3 Dynamics and energy transfer for NLS/GP 

Many of the phenomena and mechanisms present in the dynamics of NLCME (2.1) are present in the 
related and more studied nonlinear Schrodinger / Gross-Pitaevskii (NLS/GP) equation: 

i9t$ = -A$4-y(a;)$ + .gl$l2$, (3.1) 

where $ = $(a;,t) is complex- valued, x S M", i e R. 5 > corresponds to a repulsive or defocusing 
nonlinear potential, while g < Q corresponds to an attractive or focusing nonlinear potential. Before 
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embarking on a study of NLCME, we briefly review results for NLS/GP. 

To fix ideas, we take V{x) to be a smooth potential well {V < 0), decaying to zero rapidly as \x\ oo. 
NLS / GP is a Hamiltonian system with conserved Hamiltonian energy 

n[U]= [ (\VU\' + V{x)\U\' + ^-\Uf)dx. (3.2) 

Solitary standing wave solutions are solutions of the form: ^(x, t) — \l'(x)e^'^^', from which we obtain 
the nonlinear eigenvalue problem for u{x) 

= -A^ + y(a;)* + gl^-pvE-; ^ G L'^ (3.3) 

In the low intensity (linear) limit, (3.3) reduces to the the linear eigenvalue problem for the Schrodinger 
operator —A + V. 

The operator —A + V has continuous spectrum covering the non-negative half-line: > and, in 
general, a finite number of negative discrete eigenvalues fio* < ^i* < ■ ■ ■ < fin* < 0, with corresponding 
eigenfunctions Vj*) j = 0, ...,n, ||V'j*||2 = 1; recall subscript asterisks denote solutions in the linear 
limit. The cigcnfimction ^'o*(a;) corresponding to Qq* is the ground st,ate, a minimizer of the (linearized) 
energy 7Yiinear[C^] = / (| VJ/p -|- y(a;)|[/p) dx, subject to the constraint: \\U\\2 = 1. Nonlinear defect 
mode families, (^„j (a^), (a;)), bifurcate from the zero state at the linear eigenfrequencies [60]. The 
nonlinear ground state, {'^aa, ^ao) can also be characterized variationally as a minimum of 'H[U] subject 
to the constraint {[UW^ = ex.. For small L^-norm we have for j = 0, 1, . . . and a ^ 

■^^^{x)=a{i>,,{x) + 0{\g\ |a|2) ) 
n,(a)=17j. + gc%\a\'' + 0(/|a,n, a,- g C 

where c^^ > is a positive constant. 

Soffer and Weinstein [67, 68] studied the dynamics of the initial value problem for NLS/GP (3.1) in the 
case where the potential well, V{x) , supports exactly two eigenstates, "linear defect modes" , (mq* {x) , f^o* ) 
and (ui*(x), fii*). By the above discussion, there exist two bifurcating branches of nonlinear boimd states 
aji^aj)i j = 0) 1- These are used to parameterize general small amplitude solutions of NLS/GP: 

^{X,t) ~ *ao(t)(a;) + *ai(t)(a:^) + *dispersive(a;,t) 

and it is shown that generic finite energy solutions of the initial value problem converge as t ^ ±oo to a 
nonlinear ground state: ^{x,t) ^E'^.i, G C, locally in L^. This (/rownii siaie seteciion phenomenon 
has been experimentally observed; see [49]. 

The very large time dynamics, which involve energy leaving the excited state and being transferred 
to the ground state and radiation modes is governed by nonlinear master equations 

'^^"^'^ ^TPlit)Po{t) 



dt 
dPijt) 
dt 



-2TPl{t)Po{t). (3.4) 



Here, Pj{t) ~ |aj(t)p and V = 0{g'^) is non-negative and gcnerically strictly positive, provided an 
arithmetic condition implying coupling to of the discrete to continuum radiation modes at second order 
in g: 

201, - ^Q* > 0. (3.5) 

The expression for F is a nonlinear analogue of Fermi's golden rule (see also [66]), which arises in the 
calculation of the spontaneous emission rate of an atom to its ground state [17]. 

The resonance condition (3.5) also appears in a natural way in the linearized (spectral) stability 
analysis of nonlinear defect modes. In particular, consider the linearization about a nonlinear excited 
state. This yields a linear spectral problem of the form: a^HY = (3Y, where i? is a two by two self- 
adjoint matrix operator and CT3 is the standard Pauli matrix. Now H = Hq + W, where W tends to zero 
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Figure 3.1: The projections on the ground state (black) and excited state (gray, dashed) of numerical 
solutions to equation (3.1). In (a), the potential is defined with L = 1 leading to a resonance. In (b), 
the potential is defined with L = 2, leading to no resonance between the two modes. Plotted are not the 
amplitudes of projections onto the two modes, but the L^-norm of the even and odd projections of the 
solution, which are less noisy. 





quadratically in the nonlinear excited state amplitude. Since W is spatially localized, the continuous 

spectrum of (73 H is given by two scmi-infinitc intervals, the complement of an open symmetric interval 
about the origin. Now under the resonance condition (3.5), o^Hq has an embedded eigenvalue within 
the continuous spectrum. The perturbation theory of embedded eigenvalues is a fundamental problem in 
mathematical physics; see, for example, [59, 65]. These embedded eigenvalues can be shown generically 
to perturb for W arbitrarily small to complex eigenvalues, corresponding to instabilities [18]. 

In subsequent sections we explore the analogous picture for NLCME. In particular, in section 4, 
we obtain an analogous (more complicated) resonance condition, yielding embedded eigenvalues for a 
spectral problem, perturbing to instabilities and analogous dynamics of energy transfer among modes. 

We conclude this section by considering numerical simulations for NLS/GP, with a simple family of 
potentials that support two linear defect modes is 

Vl{x) = -2 sech^ {x-L)-2 sech^ {x + L). 

More detailed numerical simulations and their interpretation in light of [67, 68] is given in [62] . 

A potential well V = — 2sech^x has a single localized eigenfunction u = sechx with frequency 
17 = —1. For all values of L, the potential Vl{x) supports exactly two stationary solutions. For L 
sufficiently large, the normalized ground state is given by ■^0* ~ 2~2 (sech {x — L) + sech {x + L)), with 
frequency ilo* « — 1 — e{L) where e{L) is positive and exponentially small in L. The excited state 
is Vi* ~ 2~2 (sech {x — L) — sech (./; + Lj), with frequency fii* « — 1 + e{L) . As i is decreased, Oo* 
increases and Sli* increases toward —1. For L < 1.13, condition (3.5) is satisfied. 

We have simulated the initial- value problem for this problem with initial conditions given by a linear 
combination of the two eigenfunctions, in both the resonant and non-resonant cases. The results of these 
simulations is shown in figure 3.1 and the behavior in the two cases is strikingly different. Subfigure (a) 
shows the case L=l, in which the ground state grows slightly, while the excited state decays. In addition 
to this transfer of energy, there is a fast oscillation in the amplitudes of the two modes. Subfigure (b) 
shows the nonresonant case L ~ 2. Here there is a much larger amplitude oscillation between the 
two modes, but none of the one-way energy transfer. In section 6, we perform analogous numerical 
experiments, which confirm analytical / numerical predictions of section 5. 
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4 Linearized stability analysis of nonlinear defect modes 



We now study the linearization of the variable-coefficient NLCME (2.1) about a nonlinear defect mode. 

In the introduction, wc reviewed one type of instability that may arise. The analysis of subsection 4.2 
is focused on this scenario — instabilities arising from perturbations of eigenvalues, embedded within the 
continuous spectrum. We first discuss the conditions under which, in the limit of vanishing-amplitude 
defect modes, the linearization contains discrete eigenvalues embedded in the continuous spectrum, and 
then develop a time-dependent perturbation theory that shows when these embedded modes may lead 
to instability. In section 4.3, we discuss other types of instability that may arise. 

4.1 Conditions for embedded eigenvalues in the linecirization 

Letting E = (f^)e~*"* be a solution to the nonlinear eigenvalue problem (2.6) constrained to have 
intensity /, we linearize about E by letting 

E+ = {£+ + yi{Z, T))e-'-^ E*^ = {£*^ + ys{Z, T))e'^^ 

E_ = {£_ + y2{Z, T))e-"^^ E*_ = {£*_ + y^Z, T))e^^'^ . 

Letting yi{Z, T) = yi{Z)e~^^'^ yields the eigenvalue problem for (3 and an L"^ function y= (yi, 2/2,2/3, Vi)^ '■ 
(3y{Z) = -^^■(v{Z)+u; + ^(^^ J 5z + (^J k{Z) + W {Z)j y (4.1) 

where S3 = ( q is a 4 x 4 Pauli matrix. The matrix multiplication operator W{Z) is given by 



W 



( 2|f|2 2£+£*_ £l 2£+£^\ 

2£l£_ 2\£\^ 2£+£_ £l 

£1^ 2£X£*_ 2\£f 2£l£. 

\2£l£*_ £*J 2£+£*_ 2\£\^.) 



(4.2) 



An eigenfunction is a square-integrable solution of (4.1), with corresponding discrete eigenvalue /3. Thus 
if any non-trivial solutions exist with 3/3 > 0, the nonlinear defect mode is unstable. 

As ||^(j)|| \ 0, the spectrum of the linearization about a given defect mode approaches 

SpeCiinoarizcd = {±(t^J* - ■ ^* £ SpCCii„ear} 

In particular, from the form of equation (4.1) in the limit of \Z\ 00, the branches of continuous 
spectrum must be given by 

speci = {/3 : |/3| > Koo - \<^j*\} 

and 

spec2 = {P : \(3\ > Koo + 

so that spectrum has multiplicity two for > Kqo + In addition the discrete spectrum is given 

by 

SPeCdiscrete = {/^jtfe = - Wfe*) : 1 < fc < TV}. 

This implies that the spectral gap is the interval 

gap = (-/too + \^^j*\,Koo - |wj*|) (4.3) 
This is summarized in figure 4.1. If for some fc e {1, . . . , N}, 

\Ptk\>Woo\-M. (4.4) 

then the discrete eigenvalue P^f. is embedded in the continuous spectrum. We will refer to the boundary 
between the band gap and the multiplicity-one continuous spectrum as the "primary band edge" and the 
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Primary band edge 

Secondary band edge 
p:^. in gap Pj^;^ embedded 

J,K J,K 

ure 4.1: A schematic depicting the spectrum of the hnearization about the defect mode at zero amp 
tude. 



boundary between multiplicity-one and multiplicity-two continuous spectrum as the "secondary band 
edge." 

Example 1 For "even" defects of the form (2.9) with exactly two discrete eigenvalues W-i* = —u)u, 
therefore there will exist an embedded frequency if wi* — u)-i* = 2wi, < Kqo ~ ^i*, i-s- «/3a;i* > k^o- 
Since Koo = 'i- in this example, the condition for an embedded eigenvalue is wi* > 1/3. 

A similar calculation shows that for "odd" defects of the form (2.10), the linearization about the 
"excited states," £(±1)* always produces embedded frequencies, while the linearization about the "ground 
state" £(0)* has an embedded frequency ii oj/k > l/\/56. Under perturbation, the frequencies in the gap 
will generically remain real at small amplitude, while the embedded frequencies will be shown to move 
off into the complex plane, giving growing modes. 



4.2 Perturbation theory of embedded eigenvalues — time dependent approach 

In example 1, we observe that there are cases when the linear spectral problem (4.1), in the limit 
of vanishing nonlinear defect mode amplitude, has eigenvalues embedded in the continuous spectrum. 
The perturbation theory of embedded eigenvalues in continuous spectrum is an important fundamental 
question in mathematical physics [59, 65]. In the self-adjoint case, such embedded eigenvalues perturb 
to resonances, time-decaying states; self-adjointness precludes instability. In the present non-self-adjoint 
setting, instability cannot be precluded. 

We sketch a time-dependent approach to the theory of embedded eigenvalues [65, 67, 68] to study 
instabilities which arise. A time-independent approach can also be applied; see [18]. In particular, we 
sketch a proof of the following 

Proposition 1. Consider the linear spectral problem (4-.1), associated with a branch of nonlinear bound 
states , which bifurcates from the zero solution at frequency cuj . If the linearized operator corresponding 
to the zero amplitude limit (the operator on the right hand side of (4-1) with W = Q) has an embedded 
eigenvalue in its continuous spectrum, (see condition (4-4) )j then for arbitrary sufficiently small nonlinear 
bound states of order \\^\\, the linearized evolution equation has an exponentially growing solution, with 
exponential rate e^*, where F > (generically > 0) is of order and is given by the expression 

(4.13). 

Recall that the nonlinear coupled mode equations arc a system of scmilincar evolution equations for 
complex wave amplitudes E+ and E-. Using (2.12), the changes of variables E± = e^*®*^^)iJ± and 
ii{Z) = K(Z)e~^'® = w — inktanhkZ simplify system (2.2) to 

0^ + -3a.+ (_.«^) [f_) +A-(^„^_,^;,^*)(|-) =0, (4.5) 

The linearized evolution equation, governing the perturbation: 

cP = y = (j/+,y_,y;,t/*)e-'-^ (4.6) 
about a fixed nonlinear defect mode has the form 

idt<l> = i:3H cj) (4.7) 
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Here, H = Hq + W, with Hq and W self-adjoint operators. 



i^o = [',^ (4.8) 



'0 

where the 4x4 matrix Hq is naturally expressed in terms of 2 x 2 blocks, where 

ho = h-f' -1^5 ) (4-9) 
\-K*{Z) ojj+tdzj ^ ^ 

and the matrix W is given above in (4.2). 

Note that if ^ solves the spectral problem TisHip = Xip, then T,stjj solves the adjoint spectral problem, 
HTi^ip = Xip. U tp G L"^ we introduce, ip, the normalized adjoint eigenvector: 

= fx-^i / \ ' such that (V',^) = 1, 

where (•, •) denotes the inner product on C*-valued vector functions. 

Wc assume that there is a spectral decomposition associated with the operator Yj^H. If / € DiTisH), 
then / can be decomposed in terms of its bound state and continuous spectral parts: 

f = Pbf + Pcf, 

where 

k 

Pcf = [ (i'xj) ^xdX, (4.10) 

where ^ 

■^fe = / V. ,^ — V^' ^'^^ ( i'k , ipk) = 1 

( SsV-'fe , iPk ) ^ ' 

An analogous spectral decomposition holds for the operator Yi^^Hq. 

Consider the situation, in which Ss/fo has a simple eigenvalue Aq embedded in the continuous 

spectrum and corresponding eigenvector, i/iq. 

What are the dynamics for the perturbed evolution equation idt4> = Yi^ H (j>? 

Following the analysis of [65], we consider the initial value problem with data given by the unper- 
turbed state, Vo: 

i?Vo = AV-o, AoGa,(I]3i/o) (4.11) 

For W small (small amplitude nonlinear defect states), we seek a solution of the perturbed evolution 
equation 

idt(i> = EsiHo + W) cl>. (4.12) 
We show the existence of an exponential instability with growth proportional to e^* where 



(Es^OjV'o) 



Here ip^^^ is a generalized eigenfunction corresponding to the shifted frequency (4.16). If Aq is an 
embedded eigenvalue of Y^Hq with corresponding eigenfunction ipo having negative "Krein signature" , 
(S3^o,^o),thenr>0. 
We seek 

<Pit) = a(i)^o + Mt), {i'o,Mt)) = (4.14) 

Let denote the spectral projection (with respect to the unperturbed operator Y^Hq ) onto the part of 
the spectrum which is (i) continuous, (ii) bounded away from infinity and (iii) from thresholds (endpoints 
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of branches of continuous spectra). In [65] it is shown that the full dynamics are subordinate to the 
coupled dynamics of a{t) and <j)d{t) = Pf4'i{t), which are approximately governed by the system 

idtMt) = ^sHoMt) + a{t)P*^^Wi^o, (4.15) 

where 

. A. + fiiq*' . (4.6) 

Terms neglected in arriving at (4.15) can be treated perturbatively [65]. We consider the system (4.15) 
with initial data 

a(0) = 1, <Ad(0) = 

corresponding to the unperturbed embedded state. 

We determine the asymptotic (in time) dynamics of (4.15) by solving for cf>d{t) as a fimctional of a(t) 
and then substituting the result into the equation for a{t). To facilitate this, we first extract the rapidly 
varying time dependence of a{t) by setting 

A{f) = e*^»*a(f) (4.17) 
Then, the system (4.15) can be equivalently rewritten as 

d,Ait) = -i e^^o* ifo,WMt) ) (4 ,3) 

(SgV'cV'O ) 

idtMt) = ^sHoMt) + e-*^°*A(t) PfEgVFVo, (4.19) 
Now solving for ^d(t) with the initial condition (^d(O) = 0, gives by duHamel's principle 

<j)a{t) = -ie-''^"* / e-i{^sHo-Xoi){t-s) ^(^) pfSgVTV'o ds (4.20) 
Jo 

Substitution of (4.20) into (4.18) yields the closed nonlocal equation for A{t): 

(EaV'cV'o) dtA{t) = -(^W^Jo, ^-i{^,Ho-Xoi)it-s) ^(^^ pf^^w^^o ds ^ (4.21) 
Denote by af C ad^aHo) the spectral subset onto which to the operator P* projects. We now use the 
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spectral representation of T,sHo to compute the dominant and local contribution of (4.21). 



(S3V'o,^o) dtA{t) = 



Wipo 



= — lim ( Wtpo , / e 
£^0 V /„# 



- lim I WtpQ 



^ j^^ ( S3VA,e-*(^^^°-^°')(*-^) SsW^V'o ) V'A d\ A{s) ds ^ 

J # { e'^^"^°'^*"'^^A, VKV'o ) V-A rfA A(s) ds ^ 
gi(A-Ao)t ^ e-'(^-^°)>A, W'V'o ) V'A A{s) ds dX 



,j(A-Ao)t 



-i(A — Ao + «£) 



V^A, W^^^o Wa dA 



= — lim 



= - lim^ I W^o , j ^ 



V'A 



——t -VFVo Va ^^A ) A{t) 

i{X - Xo-ie) J J 



£^0 



= i P.V. / Va , ^WV-o ( W^V'o , V'A ) rfA 

JdT* V A - Ao / 

- TT y^ ( Va , ^(A - Ao) W^Vo ) ( W^Vo , V'A ) dA A(i) 

= % P.V. / K^A.VFV'o)!' dXA{t) - TT |(^^^,iyVo)|' AW- 

J at A — Ao 



This last line, 



(SsVo.V'o) dtA{f) =i P.V. / |(Va,VFV'o)|' dX A{t) - n | (V^^, W^Vo) |' A{t), 

J af A — Ao 



(4.22) 



(4.23) 



is precisely the conclusion of Proposition 1. The first term on the right hand side of (4.23) contributes 
an 0{W'^) = 0{\\£\\\) phase correction, while the second term determines the stability or instability of 
the dynamics. An explanation of the e— regularization in (4.22) is given in appendix D. 

Define the Krein signature of the unperturbed bound state (/:;o, to be the sign of (Es^Oj V'o)- If the 
signature is positive, then the dynamics the solution to (4.23) decay exponentially for f > and the 
dynamics are stable. If the signature is negative, then the solutions to (4.23) grow exponentially with 
0{W^) = (^(ll^ III) growth rate and the dynamics are unstable. These two cases correspond to the 
time-independent perturbation theory, which gives that positive signature embedded states perturb to 
resonances and negative signature embedded states perturb to genuine finite energy unstable eigenstates 
of asH; see [18]. 

The scenario discussed here is generic; if P = 0, then for typical small perturbations of the potentials 
K and V it will be non-zero. Since the continuum eigenmode V'Ao ^'^^ easily computed, it is diflacult 
to estimate the constant C in the formula 

r-ciiv^ollt 

Depending on the magnitude of this growth rate and the appearance of other instability mechanisms, 
discussed next, this may or may not be the most important instability of a nonlinear defect state. 
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Figure 4.2: (a) Instability scenario 2. (b) Instability scenario 3. (c) Instability scenario 4. (d) Instability 
scenario 5. 



4.3 Other scenarios for instability onset 

We here discuss four additional scenarios under which a defect mode may lose stability. Unlike scenario 
one, these transitions to instability arise at nonzero amplitude, and cannot be determined immediately 
from a condition such as (4.4) derived entirely from quantities known exactly in the zero-amplitude limit. 
In the first two scenarios, instability arises due to bifurcations involving the continuous spectrum while 
in the last two, instability arises from bifurcations involving only frequencies in the discrete spectrum. 

In a second scenario, sec figure 4.2(a), SsiJo has a symmetric pair of eigenvalues in the gap, and as 
||W^|| is increased the corresponding pair of eigenvalues of S3(iJo + W^) collide, for some critical positive 
amplitude, with the symmetrically located spectral band/gap edges. As ||W|| is further increased, each 
collision gives rise to a pair of complex conjugate eigenvalues. Those in the upper half plane correspond 
to solutions which grow exponentially as time increases, and have a growth rate given by the imaginary 
part of /?. 

In a third scenario, see figure 4.2(b) as \\W\\ is increased, a quartet of complex eigenvalues (the four 
values±/3 and ±/3*) bifurcates from the secondary band edges on both sides of the spectral gap, resulting 
in two growing modes and two decaying modes with growth rate given by the imaginary part of /3. ^ 

In a fourth scenario, see figure 4.2(c), Sa-ffo bas a symmetric pair of eigenvalues within the gap, 
and as increases, the associated discrete eigenvalues of Saif = S3(i^o + W^) collide at zero at 

some critical positive amplitude and then symmetrically ascend/descend the imaginary axis as is 
further increased. This scenario is usually associated with Hamiltonian pitchfork bifurcations in finite 
dimensions. 

A fifth scenario see figure 4.2(d) — is related to Hamiltonian Hopf bifurcations. Suppose that at 
small values of the linearized operator 'Ss{Ho + W) has two pairs of real frequencies ±/3o and ±/3i 

in the spectral gap, such that they collide pairwise at some critical amplitude. Then above this critical 
amplitude, the two pairs of real frequencies are replaced by a quartet of complex frequencies, two of 
which lead to exponential growth. As we did not investigate any defects that have more than three 
linear modes, this scenario was not observed. 



5 Numerical computation of discrete spectrum of T^^H via the 
Evans function 

We now turn to numerical computation of the eigenvalues of equation (4.1) to confirm the above analysis 
and estimate the growth rates predicted by equation (4.23). The simplest approach to numerically solving 
the eigenvalue problem (4.1) would be to discretize the derivative and to approximate the solution (4.1) 
by the solution to the large system of linear equations so derived. This requires truncating to a finite 
domain and applying artificial boundary conditions, usually Dirichlet or periodic. 



■* Closely related to this is an instability scenario described by Kapitula and Sandstede for a different system of coupled- mode 
equations. An edge bifurcation may occur when the primary and secondary band-edges coincide, i.e at the value of \\£{i) \\ for 
which (jji = 0; see equation (4.3) [39]. This was not observed in our numerical simulations. 
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Barashenkov and Zemlyanaya performed sueh a calculation using a Fourier representation in con- 
sidering the closely-related question of stability of the gap solitons of (2.7) without defect [5, 6]. This 
produced a significant number of spurious unstable eigenvalues with relatively large (0(10"^)) imaginary 
parts, and these errors were found to decay rather slowly as the number of modes used in their calculation 
was increased, possibly due to the artificial boundary conditions. Sandstede and Scheel have analyzed 
the effect that such boundary conditions has on the spectrum of linear operators [61]. Another approach 
was taken by Malomed and Tasgal, who studied the problem using averaged Lagrangian techniques to 
derive ordinary differential equations governing the evolution of perturbations to the gap solitons [48], 
although they, too, found that their method could yield spurious growing modes. Stability of nonlinear 
defect modes is studied by [46] , who report simulations of the initial- value problem, showing that certain 
nonlinear defect modes localized at a point (Dirac delta) defect are "semi-stable." 

A better way to test for instability is to use the Evans function, the analog of the characteristic 
polynomial for the linear operator (4.1). It was introduced by Evans [24] and generalized by Alexander, 
Gardner, and Jones [3]. Evans functions have since become a standard tool in both the analytic [26, 
39, 45, 56] and numerical [9, 11, 12, 13, 45, 54, 55] study of wave stability. Derks and Gottwald have 
applied the Evans function numerically to the problem of gap solitons and were able to avoid finding 
spurious unstable eigenvalues [20]. 



5.1 Definition of the Evans Function 

The Evans function is the analog of a characteristic polynomial for a class of infinite-dimensional op- 
erators. It is an analytic function whose zeroes and their multiplicity correspond to eigenvalues of 
the operator and their multiplicities. To construct the Evans function, it is convenient to rewrite the 
eigenvalue problem (4.1) as: 



i{V{Z) + (3) i {k{Z) + 2£+£l) i£% 2i£+£_ 

-i{K,{Z) + 2£l£_) -i{V{Z)+f3) -2i£+£_ -i£l 

-i£f -2i£X£*_ -i{y[Z)-ff) -i{n{Z) + 2£l£. 

2i£l£*_ i£*J i {k{Z) + 2£+£l) i{V{Z)-p) 



y (5.1) 



where V{Z) = V{Z) + uj + 2{\£+{Z)\^ + \£-{Z)\^). 

The Evans function D{(3) is constructed from the stable and unstable manifolds of the trivial solution 
y = to equation (5.1) for fixed /3 and \Z\ oo. First, rewrite (5.1) in the general form 

^y = A{Z,f3)y (5.2) 

for y G C" and A{Z, (3) an analytic n x n matrix-valued analytic function of Z and f3 which approaches 
a constant value Aoo{(3) as \Z\ oo. Each discrete eigenvalue of (4.1) corresponds to a solution of (5.2) 
that decays as \Z\ — * oo. Solutions that decay as Z ^ — oo must approach zero in a direction tangent to 
the unstable subspace of AodP) (the span of the eigenvectors whose eigenvalues have positive real part), 
which we may assume has dimension fc, while solutions that decay as Z ^ oo approacli zero tangent 
to the (n — A;)-dimensional stable subspace of {P) (corresponding to eigenvalues with negative real 
part). An eigenfunction of (5.2) exists if the stable and unstable manifolds have nontrivial intersection. 
The Evans function -D(/3) is defined as follows. Assume the matrix AooiP) has eigenvalues satisfying 

SRAi > 5RA2 > . . . > SRAfc > > KAfc+i > . . . J?A„ (5.3) 

with corresponding (generalized) eigenvectors Vi. For j = 1, . . . ,k, define Vj^i^) to be a solution to (5.2) 

satisfying the asymptotic condition rf^{Z) ^ Vj exp \jZ as Z —00. For j = k+1, . . . ,n, define rfj [Z) 
to be a solution satisfying the asymptotic condition 'n~{Z) ~ Vj expXjZ as Z ^ 00. Note that Xj, and 
r]f all depend on the frequency parameter /?. 

If a discrete eigcnmode exists with complex frequency /3, then it must approach the space spanned 
by {r]'^{Z) : j = 1, . . . ,k} as Z —00 and analogously to {t]J{Z) : j = k + 1, . . . ,n} as Z ^ 00. Thus 
these two subspaces must be linearly dependent, which will cause the Wronskian determinant 

W{Z;P) = det [ri+{Z) ... rj+{Z) n^^,{Z) ... rj-iZ)] (5.4) 
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to vanish for all Z. The eigenvectors Vj of A^^ depend analytically on /3, up to a constant multiple. We 
define the Evans function by the Wronskian evaluated at some point Zq which may be taken to be zero: 

i)(/3)=W^(0;/3). (5.5) 

The domain of the Evans function consists of the set of complex frequencies /? for which the asymptotic 
eigenvectors yy^ can be defined analytically and the condition 

is maintained. That is, the ordering of the other eigenvalues may change, so long as the (fc, n — k) 
splitting is preserved. The eigenvectors may be normalized analytically such that 

lim D(f3) = 1 

I/3HOO 

when the limit is taken along any path inside the domain of D{f3). 



5.2 Numerical experiments 

The principal way Evans functions are used numerically takes advantage of their analyticity — for a given 
closed curve 7, the winding number of £'(7) equals the number of zeros, counting multiplicity, of D{/3) 
inside 7. Since D{/]) 1 as ^ 00 in the domain of D{P), the number of zeros of D{p), counting 
multiplicity, above a given line in the half-plane is equal to the winding number of its image under 
/3^£>(/3). 

We have used winding number calculations to count eigenvalues, followed by root-finding to trace 
the evolution of these eigenvalues as the amplitude of a defect mode is increased. We discuss here some 
numerical results for defects with one, two, or three bound states in the linear regime. We looked at a 
variety of defects beyond those discussed here, with similar results, although no thorough attempt has 
been made to thoroughly explore the parameter space. 



Experiment 1: A defect supporting one bound state 

The defect defined in equation (2.10) with n = 1 supports a single linear bound state. To discuss the 

system at zero amplitude, we refer to figure 4.1. The only discrete eigenvalue is given hy (3 = luq ^ luq = 
given by a 'x' sign in the figure. Modes like those marked with the '•' and are absent for this 
particular defect. As the amplitude is increased, the eigenvalue at zero does not move. The spectral 
gap is given by formula (4.3) with wo* replaced by loq, the frequency of the nonlinear bound state. As 
seen in figure 5.1 (i), the frequency moves toward the left band edge as the amplitude is increased. In 
figure 5.2, we plot the numerically calculated eigenvalues of this defect with w = fc = 1. Figures like 
this one are repeated in this section, so it is worth spending some time dissecting it. As the defect 
mode's amplitude is increased, its frequency decreases from lj — 1 to lj — — V2 = —Kod, the left band 
edge. Note that equation (4.3) implies that the width of the gap approaches zero as loq \ —/too- The 
spectrum is symmetric across the real and imaginary axes, so only the first quadrant and its boundary 
are plotted. In the left figure, the mode's amplitude is plotted on the a:;-axis, and the real part of the 
spectrum on the y-axis. The region of multiplicity-one continuous frequencies is shaded light gray and 
multiplicity two in dark grey, with dashed lines indicating their boundaries — the band-edges — moving 
according to (4.3) and cross when ujo = 0. 

The linearization has no discrete modes of non-zero frequency at zero amplitude. Since zero remains 
an eigenvalue for all amplitudes of Sq, instabilities may only develop from edge bifurcations. Two discrete 
modes, labeled (a) and (b), bifurcate from the primary band edge. At higher amplitude, each collides 
with the band edge to produce an instability, as described by scenario 2. A third instability, marked (c) 
arises due to an edge bifurcation from the secondary band edge, in accordance with scenario 3. Mode 
(b) is the first to become unstable and develops the largest growth rate of the three modes. Below the 
large intensity 1.2, there is no instability. It is worth noting that edge bifurcations only take place on 
the band edge that moves away from the origin as ||fo|| is increased and not from the other edge. This 
same pattern is seen in all the simulations. 



19 




Figure 5.1: (i): The amplitude vs. frequency of the single nonlinear defect mode of the "odd defect" 
of equation (2.10) with co = k = n = 1. The stability of this mode is analyzed in Experiment 1 and 
figure 5.2. (ii): The "resonant", "even" defect analyzed in Experiment 2.1, with curve (a) analyzed in 
figures 5.3 and 5.4 and curve (b) analyzed in figure 5.5. (iii): The "nonresonant" , "even" defect analyzed 
in Experiment 2.2, with curve (a) analyzed in figure 5.6 and curve (b) analyzed in figure 5.7. 




ll4)ll 

Figure 5.2: Experiment 1: The eigenvalues of the linearization about the nonlinear defect mode of a one- 
mode defect, (i): The real part of the frequency, showing three modes, marked (a), (b), (c) bifurcating 
from the edge of the continuous spectrum. These lines are drawn dash-dot if the frequency is purely real 
and solid if it has a nonzero imaginary part (growth rate). The dashed lines mark the band edges and the 
light and dark shaded regions the multiplicity-one and multiplicity-two regions of continuous spectrum. 
(ii): The growth rates of the three instabilities, labeled as in (i). 
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Figure 5.3: Experiment 2.1, annotated as in figure 5.2. The frequency (i) and growth rates (ii) for the 
linearization about the left nonlinear defect mode with b = I. The embedded frequency develops a small 
instability as the amplitude is increased. A second slightly unstable mode (b) bifurcates from the band 
edge at higher amplitude. 

Experiment 2: Defects supporting two bound states 

We consider two defects in the family (2.9), with A; = 3 and different values of the parameter b. Because 
these defects support two linear bound states, the linearization has discrete consisting of the zero eigen- 
value and of another family eigenvalues that bifurcate from /? = ±2wo* which we track as ||f±i|l2 from 
zero. In experiment 2.1, this eigenvalue is embedded in the continuous spectrum like the mode marked • 
in figure 4.1. According to scenario 1, detailed in section 4.2, we expect an exponential instability for 
any ||^±i||2 > 0. In experiment 2.2, the eigenvalue is in the gap, like the mode marked + in figure 4.1. 
Surprisingly, a much stronger instability is found in the second case. 

Experiment 2.1: A defect with an embedded eigenvalue in its linearization 

In this case, we set 6 = 1 in equation (2.9) and find the linear frequencies are a;±i* = ±0.76, embedded 

in the continuous spectrum. We look first at the "left mode," for which uj-i^, — —0.76, the left branch 
figure 5.1 (ii). The real and imaginary frequencies of the linearization are given as functions of the 
amplitude in figure 5.3. As the mode's norm increases, its frequency approaches —1, so the band gap 
closes and the edges of the secondary bands go to ±2. The mode labeled (a), which is embedded in 
the continuous spectrum (like the eigenvalues marked by '•' in figure 4.1) develops an instability in the 
manner of scenario 1 (figure 1.2). The real part of the frequency does not change significantly: from 
1.52 to 1.505. A second frequency (b) bifurcates from the secondary band edge at ||5_i*||2 — 0.88 under 
scenario 3. Figure 5.3 (ii) shows that 9/3 for mode (a) is never larger than .004 and for (b) is never larger 
than 1.6 x 10~^, so no large instabilities arise. Figure 5.4 shows via a log-log plot hat the growth rate 
scales as |lf II2 for small amplitudes, in agreement with equation (4.23) (recall ||T4^|| = Odl^jp). 

Figure 5.5 shows the real and imaginary parts of the frequency for the mode with linear frequency 
Wi* = -1-0.76 (the right branch in figure 5.1). In this case, the frequency Wi decreases monotonically 
from 0.76 to -1.0 and thus the band gap widens imtil wi = and then contracts. This allows for 
some more interesting behavior. The discrete eigenvalue marked (a) begins in the continuous spectrum 
at II fill 2 = and its imaginary part grows in accordance with scenario 1. At larger amplitude it is 
subsequently absorbed in the band edge at ||5_i*||2 ~ 0.88 (essentially scenario 3) in reverse. A second 
mode, marked (b) bifurcates from the primary band edge soon after. This frequency must be real unless 
it collides with another eigenfrequency or band edge. However it soon collides again with the primary 
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Figure 5.4: Experiment 2.1: Verification that the growth rate scales as ||^||| as predicted foUowing equa- 
tion (4.23). 



band edge and picks up a nonzero imaginary part, as in scenario 2. At even higher intensity, a third 
discrete mode (c) bifurcates from the secondary band edge with nonzero imaginary part (scenario 3). 
The growth rates for this defect mode all remain relatively small, with the largest growth rate arising 
from branch (c) at larger amplitudes. 

Experiment 2.2: A defect with no embedded eigenvalue in its linearization 

For the linear defect defined by (2.9) with 6 = 3 and k = 1, the linear frequencies decrease to u)±i^ = 
±0.221 which puts the frequencies of the hncarization inside the gap. as discussed in Example 1; see 
figure 5.1 (Hi). The one discrete eigenfrequency (and its mirror image) present in the limit of zero 
amplitude is now in the gap (cf. the eigenvalues marked by '+' signs in figure 4.1), and Hamiltonian 
symmetries prevent it from developing a nonzero imaginary part without colliding with another frequency 
or band edge. We first look at the linearization about the mode with negative eigenvalue; see 

figure 5.6. In this particular example the defect mode loses stability in a collision; frequencies labeled 
(a) at /? = ±0.442 move toward zero and collide at amplitude 0.81, and, as described by scenario 4, split 
into a pair of complex conjugate pure imaginary eigenvalues, leading to an instability which reaches a 
maximum growth rate of about 0.66. This is the most significant instability seen in the three example 
defects explored. Two further instabilities, labelled (b) and (c), arise from (scenario 3) edge bifurcations 
from the secondary band edge. 

Finally, we look in figure 5.7 at the stability of the mode whose frequency in the linear limit is +0.221. 
In this case the cigcnvahie (a) in the gap moves toward the band edge, and eventually disappears in an 
edge bifurcation, moving to another "sheet" and becoming a resonance frequency with negative imaginary 
part and corresponding resonance mode (non-i^) with exponential decay in t. Another instability (b) 
arises due to an edge bifurcation when the nonlinear mode has norm about 1.25 and reaches a maximum 
growth rate of about 0.0065. 

Experiment 3: A defect supporting three bound states 

In the previous subsection, we were able to verify that the growth rate of nonlinear modes bifurcating 
from embedded frequencies scales as the fourth power of the norm, as was found analytically for 
scenario 1. In all the examples cited above the constants in front of the \\£\\'^ term was small, meaning 
that at small amplitudes, energy does not move quickly from one mode to another. Here we show an 
example where the constant of proportionality is significantly larger. 
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Figure 5.7: Experiment 2.2, annotated as in figure 5.2: The frequencies (i) and growth rates (ii) for the 
hnearization about the right nonhnear defect mode with 6 = 3. 



We consider the odd defect of the form (2.10) with oj = 1, k = 1, n = 2. The linear problem has 
discrete eigenfrequencies wq* — 1, uJi* — 2, and uj^i^ = —2. The nonlinear modes for this defect are 
described by figure 2.1 The linearizations about all three modes at zero amplitude possess embedded 
frequencies, so all three modes are unstable. We will consider the instability of the mode £i, that with 
frequency wi* = 2 in the linear limit (refer to figure 2.1), whose roots are plotted in figure 5.8. The 
mode labeled (a) bifurcates from = cji* — cj_i* — 2— (—2) = 4 moves to the other sheet in an edge 

bifurcation at i^-norm about 0.84. The mode labeled (b) bifurcates from (3^^ = oji* — CiJo* =2 — 1 = 1 
and develops a large growth rate with a maximum value of about 0.67, the one example we have found 
of an embedded mode leading to a large instabiUty. The mode labeled (c) appears in an edge bifurcation 
from the primary band edge and, as in scenario four, collides with its mirror image at the origin when 
\\£\\ « 1.43, becoming a pair of pure imaginary eigenvalues, which collide and become real again at 
amplitude 2.5 (not shown). This instability is of the same order of magnitude as branch (b). 

6 Dynamics of competing discrete modes: Time-dependent nu- 
merical simulations 

In the previous two sections, by studying the spectral problem (4.1), we provided analytical estimates 

and numerical approximations to the rate at which perturbations to given nonlinear defect modes grow. 
Since the total intensity (squared norm) of a solution is preserved, growth of perturbations can be 
expected to lead to decay of the defect mode. We expect to observe a variant of the ground state 
selection scenario of [67], outlined in section 3; energy of an imstable mode is transferred to a preferred 
or selected discrete mode and to radiation modes. This scenario for NLCME is supported by simulations 
in [29], reprinted in figure 2.1. We now present further numerical evidence for this scenario. 

At zero amplitude, the linearization about the mode has an eigenvalue /3^^ = ±(wj» — Wfe*) 

whose eigenvector corresponds to a perturbation of the mode £q^^ in the direction of the mode £(k)*\ see 
section 4. If this frequency is embedded in the continuous spectrum of ^sH -condition (4.4)-and thus 
perturbs to one with positive imaginary part, then by equation (4.23), we expect the projection of the 
solution on the mode £(^k)» to grow. Most of our simulations will begin with an initial profile described 
by a linear combination of linear defect modes. We plot, as time series, the projections of the numerical 
solution onto these modes. In section 4.2, it was shown that embedded frequencies in the linearization 
generically perturb to isolated frequencies having nonzero imaginary parts. In particular, equation (4.23) 
shows that the amplitude of the perturbation grows at the exponential rate tt | [tp^^ , Wtpo) p/| (Sa^/io, ^o) |- 



24 



0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 

ll^ll ll^ll 

Figure 5.8: Experiment 3, annotated as in figure 5.2: The frequency (left) and growth rate (right) of 
discrete modes of the hnearization about the mode £i of (2.6). 



Wc note that since the operator W scales as l|f Ip, the growth rate is proportional to \\£\\^, with an un- 
determined constant of proportionality. We have no a priori knowledge of the inner product (V'ao' ^V'o) 
(its calculation requires the continuum (non-L^) eigenfunction at frequency Aq) and the numerical cal- 
culations of section 5 demonstrate that the growth rate may in fact be quite small. 

With that in mind, we now report on the results of direct numerical simulations of the dynamics of 
the competition between defect modes. We verify the analytical and numerical predictions of sections 4 
and 5, and use the results of these sections to explain the results of [29], summarized in figure 2.2. 

The PDE is discretized using a sixth-order upwind finite-difference scheme for the spatial derivatives 
and the fourth-order Runge-Kutta method is used for time-stepping. The solutions lose energy to 
radiation through the endpoints of the computational domain, so outgoing boundary conditions are 
required, and we implement them using the method of perfectly matched layers (PML), as introduced 
by Dohnal and Hagstrom for coupled-mode equations [22]. For large- amplitude solutions, which shed a 
large quantity of radiation, the use of PML was absolutely vital to ensuring that the computations were 
not polluted by numerical artifacts. 

Experiment 1: Dynamics in a NLCME defect supporting one bound state 

In experiment 1, the defect supports only one defect mode; there is no second mode to which the solution 

may transfer energy. This is analogous to a situation studied in detail for NLS / GP [63, 64, 57, 33]. 
Our computations with the Evans function show the sole nonlinear bound state loses stability at an 
amplitude of just above 1.2 via a collision of discrete real eigenvalues with spectral gap edges (figure 5.2, 
scenario two). 

To test this prediction in a time-dependent numerical simulation, we begin with initial conditions 
consisting of a nonlinear defect mode (obtained by the methods described in appendix C) and a localized 
perturbation of about 10% of its magnitude. Wc show the results of two such experiments, the first 
with a nonlinear mode of amplitude (L'^-norm) 1.0 and the second of a nonlinear mode of amplitude 1.6 
in figure 6.1. In the first case, the system quickly sheds a small amount of radiation and then settles 
down into a stable spatial profile (with an oscillating phase not shown). In the second case, the initial 
shedding of radiation is more dramatic and a very pronounced oscillation in the shape and amplitude 
of the trapped mode is seen throughout the length of the simulation, along with a continual shedding 
of radiation. The perturbation oscillates as it grows, consistent with the quartet of complex eigenvalues 
shown in scenario two. 
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Figure 6.1: Experiment 1: In the defect supporting only one bound state, (a) a perturbed nonlinear 

defect mode with amplitude 1 is stable, while (b) a perturbed nonlinear defect mode with amplitude 1.6 
is unstable. Note the continued shedding of radiation as well as the "breathing" of the pulse shape. 



Experiment 2: Dynamics of NLCME defects supporting two bound states 

In Evans function Experiment 2. wo looked at a family of defects that support exactly two defect modes. 
In part 1, the linearization possesses an embedded frequency, while in part 2 all the frequencies of the 
linearized system are in the gap. 

Experiment 2.1: Dynamics of two bound states producing an embedded frequency 

In section 4, we show that frequencies embedded in the continuous spectrum cause nonlinear defect modes 
to lose stability, but the numerical calculations summarized in figures 5.3-5.5 show that this instability 
may be very weak. We expect to see very slow energy transfer in the dynamics of this problem. Using 
the defect defined by formula (2.9) with k = 3 and b = 1 gives embedded frequencies in the linearized 
problem. Figure 5.3 shows that when = 0.75, perturbations should grow at a rate on the order 

of 10~*. Figure 6.2a shows that the growth of the perturbation (in the f direction) is indeed quite 
small. The calculations of the spectrum of the linearization about summarized in figure 5.5 predicts 
that when = 1, a perturbation in the direction of should not grow at all (there are no 

non-zero growth rates at this point.) Figure 6.2b shows that in this case the perturbation does not 
grow, and instead decays to zero rather quickly. In figure 6.2c, the two modes are initialized with equal 
nonzero amplitude, and both lose energy, with neither appearing to grow at the other's expense. 

Experiment 2.2: Dynamics of two bound states with no embedded frequency 

Here we verify the scenario four instability described in Experiment 2.2. For that defect , which supports 
two bound states, neither of which has an embedded frequency in its linearization. In the linearization 
about the bound state with positive frequency, figure 5.7 shows that no instabilities exist with sig- 
nificant growth rates. Our numerical simulations we have confirmed this. The results of these simulations 
are qualitatively very similar to those provided in figure 6.2. 

For the defect mode with negative frequency described by figure 5.6, all frequencies in the 

linearization are real and in the gap for amplitudes below about 0.8. At this point, two real frequencies 
collide at the origin and develop nonzero growth rates (scenario four). Figure 6.3 shows the evolution 
of the initial condition 0.6 times the normalized linear eigenstate (below instability threshhold) and 0.9 
and 1.0 times the normalized linear eigenstate (above the threshhold). It is clear that energy flows from 
to 5(1) only when 5(_i) is unstable as predicted by experiment 2.2 (figure 5.3) and that the rate 
of energy transfer increases as the amplitude is further increased. 
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Figure 6.2: Dynamics experiment 2.1: (a) energy initially mainly in the mode ^(-i) (case of figure 5.3). 
(b) : energy mainly in mode ^(+i) (figure 5.5); (c): equal energy initially in both modes simulated for 
a longer time. Key: in each figure, the solid line is the amplitude of £(+i) and the dashed line is the 
amplitude of 
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Figure 6.3: Dynamical simulations illustrating the bifurcation shown in figure 5.6. (a): a stable solution 
(amplitude .6) below the merger of two real frequencies creates a pair of pure imaginary frequencies and 
instability, (b) and (c): the evolution of an initial condition with L^-norm 0.9 and 1.0, respectively, above 
the instability threshhold. Here solid lines give the amplitude of £(-i) and dashed lines give £(+i)- The 
insets present the same information with a logarithmic scale on the y-axis, demonstrating absolutely no 
growth of the second mode in (a) and exponential growth in (b) and (c). 
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Figure 6.4: (a): Initially, almost all energy is in the mode £^iy In the left image, a small perturbation is 
added in the direction of mode £[-i) with little transfer of energy, (h): A perturbation is added in the 
direction of mode £^(o)) which grows in time. This is predicted by figure 5.8. 



Experiment 3: A defect supporting three bound states 

The spectrum of the linearization about the mode £(i) of a defect supporting three bound states was 
discussed in figure 5.8. This mode has frequency uj\ = 2. The other two frequencies of the linear problem 
arc luq ^ 1 and w_i = — 2. Our arguments show that the eigenmode of the hnearization with frequency 
= 1 = LOu — ciJo* corresponds to perturbations in the direction of £(o) , while those with frequency 
/3i_i — 4: — wi* — W-i* correspond to perturbations in the direction of mode £(-i)- Figure 5.8 shows 
that the perturbation in the direction of f (o), labeled (b) in the figure grows at a much faster rate than 
the pertubation in the direction of ^(-i), labeled (a). In this last time-dependent experiment, we aim 
to test this prediction numerically. 

We have performed two simulations. In the first, figure 6.4a, a perturbation of the mode in the 
£(-1) direction quickly decays, while in the second simulation, figure 6.4b, a perturbation of the mode 
ill the 5(0) direction grows, eventually overtaking the magnitude of the £(^ij mode. This is consistent 
with scenario one, as well as with the relative sizes of the growth rates found numerically in Experiment 
3. 



7 Summary and discussion 

Wave propagation through a periodic one- dimensional structure at Bragg resonance is governed by 
NLCME. The introduction of a defect leads to a generalization of NLCME having spatially varying 
coefficients (potentials), which are determined by the defect. The linear coupled mode equations with 
defect potentials may have localized eigenstates. For a nonlinear medium, these linear states persist 
as nonlinear defect modes, which bifurcate from the zero state and the linear (defect) eigenfrequcncies. 
Such nonlinear defect modes correspond to optical pulse states, pinned at the defect sites. Such periodic 
structures and the nonlinear defect states are candidates for optical storage schemes [29, 30]. Only stable 
nonlinear defect states could play this role. 

Therefore, we have considered by analytical methods and numerical simulation the stability and 
instability of such nonlinear defect modes. We first classify the various scenarios for onset of instability 
and give examples which arise for a large class of structures. 

We investigate in particular detail the situation of a family of bifurcating states, whose linearization 
has an embedded eigenvalue in the continuous spectrum, in the limit of zero amplitude, i.e. at the 
bifurcation point. Here we find (generically) that for arbitrarily small amplitudes such states are expo- 
nentially unstable with a exponential rate of growth of order the fourth power of the nonlinear bound 
state. Although exponentially unstable, we find this growth rate in practice to be quite slow in many 
examples. 
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States for which the hnoarization, at the bifurcation point, exhibits no embedded eigenvahies appear 
to be stable, at least on long time intervals, for a range of amplitudes but exhibit instability via a collision 
of discrete eigenvalues of the linearization or the collision of a discrete eigenvalue with the endpoint of 
the continuous spectrum. For a given defect with multiple nonlinear bound states, different states may 
follow different scenarios and instability may result in the transfer of energy from one mode to another 
and to radiation. 

It is instructive to compare the dynamics we have found with those of the NLS/GP equation. In the 
small-amplitude limit, the behavior is slightly different because NLCME's second branch of continuous 
spectrum makes it harder to find defects which will produce linearizations without a resonance. At 
higher amplitudes, we have found a large number of instabilities that arise due to collisions of discrete 
eigenvalues with the endpoint of the continuous spectrum. Instabilities due to collisions between pairs 
of eigenvalues lead to the strongest instabilities (largest growth rates) observed. A related problem, 
also with two branches of continuous spectrum, is a coupled NLS system describing light propagation in 
birefringent optical fibers, studied by Li and Promislow in [45]. Using somewhat different methods than 
those employed here, they also show a solitary wave is unstable due to bifurcations involving embedded 
eigenvahies. 

Another interesting behavior for the NLS/GP system has been seen by Kirr et al. [41] and some of the 
numerical experiments discussed above suggest the same phenomenon exists in NLCME. They consider 
NLS/GP system with a defect equal to the superposition of two potentials placed a fixed distance apart 
and discover a symmetry-breaking bifurcation. At small amplitudes, the stable ground state solution is 
given to leading order by the sum of two ground states of the individual potentials, in-phase and of equal 
amplitude, a highly symmetric configuration. At higher amplitudes, this configuration loses stability 
and the stable solution is concentrated almost entirely at one or the other defect location. This is what 
is known in dynamics as a supercritical pitchfork bifurcation. Bifurcation scenario four is precisely the 
mechanism seen in a generic pit(;hfork bifurcation, so it is of interest to numerically locate the stable 
states that should exist at amplitudes above the critical amplitudes at which bifurcations are seen in 
experiments 2.2 and 3. 

The symmetry breaking bifurcation could be understood by appeal to a variational principle. In 
cases where the ground state can be characterized as a minimum of Ti-lU], the NLS-GP Hamiltonian 
energy, subject to fixed ||C/||| = Af, it can be shown [4] that for small Af the ground state is symmetric, 
with peaks situated over the local minima of each well, while for A/" sufficiently large the ground state is 
strongly localized in one well or the other. Symmetry breaking, if it occurs in NLCME, would be quite 
interesting and could be studied by the methods of [41]. Note however that in contrast to NLS-GP, 
nonlinear bound states of NLCME are critical points of infinite index, since the linearized energy has 
continuous spectrum which is unbounded above and below. 

Finally, it would be instructive to extend (4.23) — which describes the growth or decay of a single 
defect mode due to its interaction with the continuum -to a coupled system of equations describing the 
interactions of multiple defect modes with each other and with the continuum. Such a derivation has 
been carried out rigorously by Soffer and Weinstein for the NLS/GP system with two defect modes. 
They show that in the NLS/GP at low amplitude, the ground state is nonlinearly stable and derive 
a Fermi golden rule that predicts that half the energy of the excited state will be transferred to the 
ground state as t oo. Immediate extensions of this are to NLS/GP with 3 or more defect modes 
[70] and to NLCME, especially in the case of 3 or more modes. Such models were able to predict the 
symmetry-breaking bifurcation in [41] and may be able to illuminate the Hamiltonian bifurcation of our 
experiment 2.2 in section 5. We note that Aceves and Dohnal have derived finite-dimensional models 
in [1], although these models do not incorporate the coupling to the continuous spectrum seen to be 
important in section 4. 
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A List of Symbols 

Complex variables: and ?sz denote the real and imaginary parts of a complex number z. Its 

complex conjugate is z* . 

A projection operator: - projection onto the continous spectral part of the self-adjoint operator H. 
Inner Product Convention: {f,g) = J f{Z) g{Z) dZ 
Pauli matrices: 

^^^=(5 o);^3=(o -l)'^3=(o ?)' where7=(^J J 

Linear and nonlinear solutions: £ = (|+) is a solution to the nonlinear eigenvalue problem (2.6) 

while solutions with subscript atcrisks, i.e. are solutions to the linear problem (2.8). Associated to 
these are frequencies oj (eigenfrequencies w* of Ssi? in the linear limit.) 
Plemelj identities 

{xTiOy^ = P.V. a;-^ ± iw S{x) (A.l) 



B Sketch of derivation of NLCME with potentials 

Consider light propagation in a fiber with Uncar polarization and assume that the signal propagates with 
a constant (given) transverse profile. The dielectric susceptibility is assumed to possess Kerr (cubic) 
nonlinearity. Under these assumptions, Maxwell's equations for the scalar electric field E{z, t) simplify 
to 

d^_E{z, t) - c-'^n^{z, E)dfE{z, t) = 0. (B.l) 

If the dielectric function is assumed to be constant n{z, E) = n, this is the linear wave equation. 
Seeking plane waves of the form e'('^^~'^*) , leads to the dispersion relation 

^2 _ ^2^2/-2 ^ Q (g_2) 

For any pair (fc, w) satisfying (B.2), the solution is a Fourier sum of backward and forward propagating 
plane wave solutions: 

E{z,t) = E+e'^^'-'^*"' + E_ e-'('=^+'^*) + c.c. (B.3) 

where E+ and E^ are the complex amplitudes of the forward and backward moving Fourier components 
at wavenumber k. 

Nonlinearity and periodic longitudinal variations in the linear dielectric constant will result in cou- 
pling between the forward and backward plane waves, and slow changes to the periodic profile will result 
in the creation of a potential. The dielectric constant consists of three parts = + «grating(^) + '^nl 
where 



'"grating 



fi^ =n'^ + eWiez); 

(z) ^2ex'-^Hez)cos{2koz + 2(j){ez)); (B.4) 
n'^^=X^'^E\z,t). 



The dimensionless parameter e is taken to be small and the remaining parameters to be order one. The 
form of the dielectric constant B.4 implies a specific relation between (a) the amplitude of refractive 
index variations, (b) the spatial scale over which the uniform grating is modulated (by apodization, 
introduction of defects etc.), and (c) the nonlinearity. This relation is chosen so that all these effects 
enter at the same order in the multiple scale analysis [40], i.e. the maximal balance. This yields an 
ansatz of the form: 

E = v^E^o + 0(e^/^) where 

Eo = £+(^i,ti)e'('=°^-'^°*+*(^i)) -|-£^_(^;i,ii)e-'(*''^+'^°*+*(^i». 
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Here z\ = ez and t\ = et. Notice the approximate solution is chosen with twice ho wavelength of the 
underlying grating — this is the Bragg resonance condition. Inserting this ansatz into (B.l), separating 
terms by magnitude in e, and eliminating secular terms at order e^/^ in the expansion leads to "slow 
equations" for the evolution of the amplitudes E+ and E- : 



c 



-id^^E_+i-dt,E_ + 



2n 



2 -<^'(-^i) 



E^ 



E_ + 



(3) 



2nc 

/I) 



E^ 



'^oX '{zi) f. 



2nc 



E+ + 



-{2\E_f + \E+\^)E+=0; 



{2\E+f + \E_f)E_ = 0. 



This multiple scale procedure [40] is rigorously justified for a close analog of these equations in [32]; see 
also Martel [50, 51] for higher-order corrections that can effect the stability. 

Introduce dimensionless distance Z = zi/L, time T = nZ/c, and electric field E± = E±/E variables 
for some characteristic length scale L and electric field strength E and scaled functions 



Wl{Z) = W{LZ), MZ) = cP{LZ), and x'iHz) = x^^KLZ). 



(B.6) 



A sketch of the above scalings and the Bragg resonance condition is given in figure 1.1; for a more 
detailed discussion of the physical length and time scales see [29]. The above equations become 

idTE+ + idzE+ + k{Z)E_ + V{Z)E+ + (jE+j^ + 2\E_\'^)E+ = 
iOtE. - idzE_ + k{Z)E+ + V{Z)E_ + {\E_\'^ + 2\E+\'^)E_ = 0, 



where we have defined 



k{Z) 



i^qL (1) 
2nc' 



which is precisely 2.1. Note that the coefficient k{Z) arises entirely due to modulations in the depth of 
the grating, while V{Z) may arise entirely due to slow changes to the average susceptibility or else due 
to local modulations in the wavelength of the grating. 



C Numerical Computation of Nonlinear Defect Modes 

Equation (2.6) is discrctized on the Z-interval [— L,L], setting h = ^ and Zj = —L + j h, for j = 
0, . . . ,N — 1. Approximating E±{Zj) = e±j, and defining the A/'-dimensional vectors of unknowns e+ 
and e_, we arrive at the system 

ive+j + i{Ve+)j + K{Zj)e-j + V{Zj)e+j + (|e+,,f + 2|e_,,f )e+j = (C.la) 
Lve-j - i{Ve-)j + K{Zj)e+j + V{Zj)e-j + (|e_,,f + 2|e+,,f )e_,,- = 0, (C.lb) 

AT 

h-J2{\e+,,\' + \e.,jf)=I. (C.lc) 
j=i 

with 2A'' complex unknowns e±j and one real unknown oj. V is the discretized derivative operator, 
which will be discussed further below. Below we use symmetries to reduce the number of real unknowns 
from 4Ar+ 1 to 7V + 1. 

Assuming k{Z) and V{Z) are even functions, the linear standing-wave solutions to (2.8) all satisfy 

where s = ±1 and the overbar denotes complex conjugation. The nonlinear defect modes are the 
continuations of the linear defect modes and preserve these symmetries, although it is possible that 
these solutions lose stability in symmetry-breaking bifurcations at higher amplitudes.. Incorporating 
this into the scheme (C.l) eliminates 2^^ complex unknowns. 
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By phase invariance, -B+(0) may be assumed to be real and positive. Then ^{E^) is even and 
is odd. Using this symmetry, we need solve only for Rj = 3?(e+j) for j = 0, . . ■N/2, and Ij = 9(e+j) 
for j = 1, . . . N/2 — 1. The remaining values may be found by symmetry, leaving N + 1 real unknowns. 
Accounting for the symmetries eliminates a one-dimensional invariant set of fixed points and speeds the 
root-finder's convergence. 

The eigenvalue problem is posed on an infinite domain. Solution of a problem on an infinite domain 
involves balancing two competing sources of error: truncation and discretization. To eliminate error due 
to the truncating the computation to an interval [—L,L], L, and thus dZ should be taken large, but 
to eliminate discretization error dZ should be taken small. To balance these two effects, a non-uniform 
discretization is imposed. The finite interval [— tt, tt) is mapped to a much larger interval by 

C = sinhAZ, (C.2) 

as proposed by Weidemann and Cloot [72] and summarized by Boyd [8] . The uniform grid Zj = —tt+JAZ 

is then mapped to a nonuniform grid with points that are closely and almost uniformly spaced near 
C = 0, and sparsely spaced for large (. Boyd provides guidelines for the choice of the stretching parameter 
A [7]. Because the Evans function requires an integration from conditions at ±oo, it is important to 
resolve the tails well at this stage. We found we could obtain more accurate numerical defect modes 
using 256 points on a stretched grid than using 2048 points on a uniform grid, where both grids were of 
the same width, in addition to the obvious savings in computation time. 

On the periodic domain [— 7r,7r), Fourier collocation is used to derive the discretized derivative op- 
erator V in (C.l), which is then multiplied by a weighting factor Wj = ^^^_^^2 ^'^ account for the 

stretching transformation (C.2). As long as the solutions are sufficiently small at Cmax = sinhATr, the 
error introduced into the equations by the assumption of periodicity is small, and the numerical solutions 
converge quickly as n, the number of grid points, is increased. What is important is that the exponential 
decay rate in the tails is computed accurately and that the numerical solution is computed accurately 
far enough into the tail so that \E±\ 10~^, then since the linearized equations will have coefficients 
depending on \E±\'^ ^ 10~^^, any error due to the truncation will be on the order of that due to using 
finite-precision arithmetic. 

System (C.l) was solved numerically using programs from the MINPACK suite [53]. For the smallest 
chosen value of J, the initial guess for the nonlinear solver was a suitably scaled linear defect mode. For 
larger values of /, the initial guess was extrapolated from previously computed solutions. The program 
continued producing new iterates until it reached the left edge of the band gap, and was able to find 
reasonable solutions quite close to the band edge. 

D Technical details from section 4.1 

We explain the manner in which the limit £ — > is taken in equation (4.22). 

^ ( e-'^^-'^°>ipx, Wipo ) V'A dX A{s) ds 
= lim / e'(^-^°+'^)* Va dX f ds A{s) ( e-*(^-^°+'^)>A, Wipo ) 

Ja* Jo ^ ' 

= iA{t) [ i dX ( ^x,W^Po ) - iAo [ Va dX ( ^a.^^V-o ) 

Ja* X- Xo+lO Jaf X — Xo+lQ 

+ * / / ^ ? ^ dX i^x, W^o ) dsA{s) ds (D.l) 

Jo Jaf X — Xo + lO 

The first term on the right hand side of (D.l) contributes to (4.23). We now show that the remaining 
terms decay. In particular, it suffices to show that for any smooth and decaying function, /, the term 

Ja* X- Xq+1\) 
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decays as t 1 +00. Let x('^) be a smooth function, which is equal to one on af, vanishes outside a 
neighborhood of af and vanishes near thresholds and infinity. For any k> 1 

r pi{X-Xo)t r gi(A-Ao+ie)t 

/* ^ X ^-n ^^ i^xJ) dX= hm / ^ r ^. (V-A ,/) i^x dX 
J a* X- Xo + lO J a* X- Xo+ie 

lira -i [ I ( Va , / ) e'(^-^°+'^)* Va dX ds 
lim -if / ( Va , / ) V'A X(A) e'(^-^''+'^> dA 



= lim ~iJ^"J{i^x,f)i^x X(A) ( (it)-'^^)' e^(^-^°+'^)« rfA 
= lim -^(-^^)-'=y" J d'^[ { , f ) i'x x(A)] e*(^-^o+fe). 

(D.2) 

E Numerical computations with Evans Functions 

Before describing the numerical method used to compute the Evans function, we first comment on the 
sources of error in this calculation. The largest source of error, is a loss-of-significance in the calculation 
of the Wronskian (5.4). As each column in this determinant is the image of a solution that grows 
exponentially as Z moves from ±00, the norm of each column vector may be quite large. Near zeros 
of the Evans function D{Z), there must be largo cancellations that arise in evaluating the determinant 
that defines it. To mitigate this loss-of-significance, the solutions to the differential equations (5.2) must 
be extremely accurate. First that means that the defect modes around which we linearized must be very 
accurate, and must bo well-computed deep into the exponentially decaying tails, to allow the shooting 
to start in a region where the linearized problem is very nearly the linearized differential equation at 
Z = ±co. Additional considerations that arise are more technical and are discussed below. 



E.l Numerical Computation on Exterior Product Spaces 

The numerical calculation the Evans function is extremely delicate, and there exists a wide literature 
discussing how to accurately implement the calculation [10, 13, 69]. We will summarize briefly and 
point to the appropriate references. The most straightforward calculation would be to first find the n 
vectors r]^ (0) andn ri~ (0) and calculate their Wronskian (5.4) directly. That is, for each of the k positive 
eigenvalues, choose L large enough that A{±L) k. and numerically integrate 

"W-"''^ (E.1) 

from Z = —L to = for j = 1 ... A; and 



^=^(^H- (E.2) 
rij \z=L = Vjt^'^ 

from Z = L (backwards) to Z = for j = /c -|- 1 . . . n. 

When the stable and unstable subspaces of are one-dimensional, this strategy works well, but 
when they are of higher dimension, this works poorly. Consider the equation 

%=AiZ)y (E.3) 
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y € C", such that the two largest eigenvalues of A arc positive with /i2 > /^i > and corresponding 
eigenvectors y2 and yi. In attempting to compute the solution y~ yie'^^'^, rounding errors will inevitably 
introduce a small error in the direction of y2 which will grow at the faster exponential rate /Lt2, which, 
for the necessarily large values of L, may cause significant errors. One might hope to correct this error 
by performing a re-orthogonalization at each step, but this, too, is problematic, as discussed in [13, 20]. 
An algorithm that overcomes these difficulties without resorting to a compound matrix formalism has 
recently been constructed by Humpherys and Zumbrun [37]. This is important for computations in C" 
with 1 since the dimension of the exterior product space grows very rapidly with n. 

A method was proposed by Brin [13] and significantly refined and made mathematically precise by 
Bridges and collaborators [10, 11, 12], working in various combinations. We briefly introduce the basic 
features of the method, and refer the reader to the cited papers. Let Sj be a set of orthonormal basis 
vectors for C™, then the exterior product (the wedge product) is defined on the Cj by the properties of 
bilinear ity: 

Bi A {asj + Bk) = aci A Sj -Veihe}^ 
{aei + ej) Aek = ae, A ej + ej A e/j 

and antisymmetry: 

which implies that A e j = 0. The exterior product is extended to C" using (bi)linearity and antisym- 
metry: if u = Y^^=i ^j^j ^^cl V = Y^^=i '^j'^h then 

ViUj)ei A Cj. 

One may easily see that the set {cj A ej} forms the basis for the vector space of 2-forms over 

C™. The 2-form Vi A Vj gives the complex area of the parallelogram with sides Vi and Vj. Given an 
i-form and a j-form with i + j < m, we may define their wedge product as an {i + j) form using the 
anticommutativity and bilinearity. This is a binary operation of the form 

A : CC^) X cC?) c(-+^). 

Note now that the Evans function can be written in exterior product notation as 

£>(/3) = T7+(0) A . . . A 77+ (0) A 77^+1(0) A ... A 77" (0) 
= >V+(0) AW-(O), 

where 'W^{Z) is the fc-form given by the wedge product of the ri^{Z) and 'W~{Z) is the [n — A:)-form 
given by the wedge product of the riJ{Z). In the present case, : 1 < j < 4 is the standard basis of 
C*, and the space of 2-forms A^(C*) is six-dimensional with orthonormal basis vectors: 

/i = 6*1 A 6*2, /2 = ei A 63, /s = ei A 64, /I = 6*2 A 63, = 6*2 A 64, /e = 63 A 64. (E.5) 

The essential idea of all the so-called compound matrix methods is to numerically evolve W^{Z) 
from Z = =fL to Z = 0, rather than integrating the individual solutions whose wedge product forms 
W+. If y^ through y*' satisfy evolution equation (E.3) then F = A . . . Ay* satisfies the linear evolution 
equation: 

k 

j=l 

fc 

= ^ jT^ A . . . A j?<^-i) A {A{Z)y>) A y^^+^^ A ... A y^. 



uAv 



n-l 

E 



E 



{UtVj 



i=l 
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Using basis vectors for the wedge space a'^(C"'), this equation may be rewritten in matrix form as 



dY 
dZ 



(E.6) 



where ^^'^^ (Z) is a (^) x (^) matrix. In the present situation, n = 4 and k = 2, then this may be written 
out, using basis (E.5): 



I an + 022 023 024 -ai3 

^32 On + 033 034 O12 

042 ^43 ail + 044 

— 031 021 022 + 033 

-O41 021 043 





-041 



0131 



-a42 



-Oi4 \ 

— Oi4 

ai2 ai3 

034 — 024 

022 + 044 023 

032 033 + 0144 / 



(E.7) 



If /Ltj and /Zj are eigenvalues of A with eigenvectors vi and Vj . Then Vi A Vj is an eigenvector of AS^^ with 
eigenvalue \i = Thus only the solution with fastest growth need be (stably) computed. 



E.2 Integration scheme 

A vector v € a'^(C") is called decomposable if it can be written as the single wedge product of fc elements 
of C", and not all such vectors are decomposable. A vector v = v^^j in A^(C'') is decomposable 

if and only if the 4- form v l\v = which gives the algebraic condition 

v\v^ — V2V5 + V3V4 = 0. (E.8) 

This quantity is the Grassmannian invariant, and the set on which it vanishes is called the Grassmannian 
manifold. The fc-form W"*" is by its definition decomposable, and it is important for numerical accuracy 
that the computed approximation to be so as well. Therefore wc require a numerical integration routine 
that conserves the Grassmannian invariant. One such family of methods is the fully implicit Runge- 
Kutta schemes using the Gauss-Legendre points [12, 38]. We use the three-step sixth-order scheme of 
this type: 

3 

r+^ = r+^zj2bjKi (E.9) 

1=1 

where AZ is the stepsize, and the superscript n denotes points in the spatial discretization. The 
increments Ki are implicitly defined by 

3 

Ki = A{Zn + CiAZ) -ir+Yl ""ii^o) (E-10) 
and the coefficients are given by: 




Since the matrix a has no nonzero elements, all the substeps Ki of the Runge-Kutta algorithm must be 
calculated simultaneously by solving (E.IO). Thus, although equation (5.1), used to construct the Evans 
function has four dependent variables, the number of unknowns at each step is increased to six by the 
wedge product formulation, and then tripled to eighteen by the numerical scheme. Since equation (E.6) 
is linear, this is still a small amount of work at each step. 
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E.3 Additional digit-preserving numerical measures 

For large values of the matrix in equation (5.1), and thus in (E.7) is nearly constant, the solution 
is not far from the product of the exponential of this matrix and the initial vector, which is especially 
simple given that the initial condition is an eigenvalue of A^"^^ . For most values of p we are interested in, 
we find that the contribution of the constant part of A dominates the computed solution so it is useful 
to solve this part exactly-the growth or rapid oscillation due to the largest eigenvalue. 

If V is an eigenvector of a matrix M, then it is also an eigenvector of M — cl for any c, so we modify 
equation (E.l) to use 

instead of A. Thus the compound matrix A^'^^ is replaced in equation (E.6) by A^"^^ = AA A = A^'^^ — 
Ai — A2, and the Runge-Kutta algorithm is only used to compute the non-constant part of the evolution. 
The solution ai Z = can then be multiplied by e(^i+^^)^ to find the solution. While the integration 
scheme converges at sixth order in AZ with or without this modification, we found empirically that the 
error found using the modified method to be much smaller for each fixed value of AZ, even for values 
of P near zero. 

An additional source of error comes from solving (E.6) from ±L to 0, rather than from ±00. We take 

L large by working on the same nommiform grid as was used in the solution to (C.l). This is especially 
important if the nonlinear mode is only known numerically. Another solution to the tail problem is to 
use a rapidly convergent series form of £*(/?) to obtain the solution of equations (E.l) and (E.2) at some 
finite value z = zkzo, and then solve the odes numerically from zLzq to [45]. 

The code was tested by confirming the results of Berks and Gottwald [20] on the stability of gap 
solitons in NLCME without defect. We found that using the non- uniform grid for integrating the 
ordinary differential equation (E.6) in the shooting method greatly improved the accuracy of -D(/3), 
especially in a neighborhood of /? = 0, where a zero of high multiplicity decreases the accuracy of the 
computation. In figure E.l, we show that the computation on the nonuniform grid with 512 points gives 
a more accurate count of the zeros than a uniform grid calculation with 2048 points which predicts two 
spurious unstable eigenvalues. The tests also showed that the interval of computation [— L, L] must be 
taken quite large to properly resolve the tail's effect on the computed D{(3). Our numerical solutions 
of (C.l) were not able to resolve the solution for Z this large without using the stretched grid. The value 
of L used in the calculations section 5.2 had to be taken fairly large before the computation converged. 
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